source("read.data.R")
source("util.R")
source("analy.model.R")
source("analy.model.selection.R")
source("analy.QRE.general.R")
source("analy.plot.R")

library(lme4)
library(lmerTest)

# e0<-par.setting[1] 
# ea<-par.setting[2]
# inv.cost<-par.setting[3]
# p0<-par.setting[4]
# p1<-par.setting[5]
# penalty<-par.setting[6]
# 
# gamma.pay<-par.behavior[1]
# gamma.r<-par.behavior[2]
# gamma.a<-par.behavior[3]
# gamma.i<-par.behavior[4]

time.start<-proc.time()

library(parallel)
cluster<-makeCluster(6)
clusterExport(cl=cluster,ls(),envir=.GlobalEnv)

cn2<-c("gamma pay","gamma ransom","gamma attack","gamma invest","follow me inv","pay social norm","pay fairness","inv social norm","follow you inv")
cn.list<-list(cn2)

treat<-"Base"
data<-read.file.index(1)
data<-subset(data,data[,"Period"]>1)
data.base<-subset(data,data[,"Treatment"]==treat) 

par.setting<-c(data.base[1,"Data1"],data.base[1,"Outside"],data.base[1,"Investment1"],0.8,0.5,data.base[1,"Penalty"],0,0,0)
par.b<-c()

par.b<-convert.par.stl(ms[[2]]$par,model.spec=ms[[1]])
par.b<-as.numeric(par.b)
par.b<-ifelse(is.na(par.b),ms[[7]],par.b)

par.b[5]<-0
par.b[9]<-0

par.b<-c(0.015,0.15,0.05,0.1,0,0,0,0,0)

# 
# r<-pmax(1,data.base[,"ransom"])
# inv<-ifelse(data.base[,"attack.decision"]==1,data.base[,"i_1"],data.base[,"i_2"])
# pay<-ifelse(data.base[,"attack.decision"]==1,data.base[,"p_1"],data.base[,"p_2"])
# 
# pp1<-prob.pay(c(1:100),1,par.setting=par.setting,par.behavior=par.b)
# pp0<-prob.pay(c(1:100),0,par.setting=par.setting,par.behavior=par.b)
# #pp.model<-mean(na.omit(pp[r]))
# #pp.fit<-mean(na.omit(c(data.base[,"p_1"],data.base[,"p_2"])))
# 
# data.pay<-cbind(r,inv,pay)
# data.pay<-subset(data.pay,(!is.na(data.pay[,"r"]))&(!is.na(data.pay[,"pay"])))
# 
# pp.data<-c(mean(data.pay[data.pay[,"inv"]==0,"pay"]),mean(data.pay[data.pay[,"inv"]==1,"pay"]))
# pp.model<-c(mean(pp0[data.pay[,"r"]][data.pay[,"inv"]==0]),mean(pp1[data.pay[,"r"]][data.pay[,"inv"]==1]))
# 
# pp.summary<-cbind(c(0,1),pp.model,pp.data)
# dimnames(pp.summary)[[2]]<-c("inv","model","data")
# 

qre.r.calc<-function(par,par.behavior=NULL) {
  
  par.b<-par.behavior
  par.setting<-par
  
  u.m1<-matrix(nrow=2,ncol=2)
  u.m2<-matrix(nrow=2,ncol=2)
  
  for(i in 0:1) for(j in 0:1) {
    u.m1[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
    u.m2[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
  }
  
  dimnames(u.m1)<-list(c("0","1"),c("0","1"))
  dimnames(u.m2)<-list(c("0","1"),c("0","1"))
  
  qre.i<-find.qre(c(par.b[4],par.b[4]),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  # qre.i<-find.qre(c(0.01,0.01),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  
  qi1<-qre.i[[1]]
  qi2<-qre.i[[2]]
  inv.summary<-c("no-no","no-yes","yes-no","yes-yes")
  inv.summary<-cbind(inv.summary,c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2]))
  inv.p<-c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2])
  
  
  qrr1<-qre.r(1,par.setting=par.setting,par.behavior=par.b)
  qrr0<-qre.r(0,par.setting=par.setting,par.behavior=par.b)
  
  qra00<-qre.att(c(0,0),par.setting=par.setting,par.behavior=par.b)
  qra01<-qre.att(c(0,1),par.setting=par.setting,par.behavior=par.b)
  qra10<-qre.att(c(1,0),par.setting=par.setting,par.behavior=par.b)
  qra11<-qre.att(c(1,1),par.setting=par.setting,par.behavior=par.b)
  
  prob.r0<-inv.p[1]*sum(qra00[2:3,2])+inv.p[2]*qra01[2,2]+inv.p[3]*qra10[3,2]
  prob.r1<-inv.p[2]*qra01[3,2]+inv.p[3]*qra10[2,2]+inv.p[4]*sum(qra11[2:3,2])
  
  prob.sum<-prob.r0+prob.r1
  
  prob.r0<-prob.r0/prob.sum
  prob.r1<-prob.r1/prob.sum
  
  er<-sum(prob.r0*qrr0[,1]*qrr0[,2]+prob.r1*qrr1[,1]*qrr1[,2])
  er   
} 

par.setting.list<-list(par.setting)
data.qre.r.penalty<-calculate.multiple(qre.r.calc,par.setting.list,key.index=6,range=c(0,100),n=20,par.behavior=par.b)
data.qre.r.penalty<-cbind(data.qre.r.penalty,100-data.qre.r.penalty[,1])
# data.qre.r.pay<-calculate.multiple(qre.r.calc,par.b.list,key.index=6,range=c(0,50),n=20,par.setting=par.setting)

# plot(data.qre.r.penalty[,1],data.qre.r.penalty[,2],type="l",xlab="penalty",ylab="expected ransom",ylim=c(0,100))

#jpeg(file="output/sim.r.penalty.jpg")
plot.multiple(data.qre.r.penalty[,1],cbind(data.qre.r.penalty[,2],data.qre.r.penalty[,3]),c("boundedly rational attacker","rational attacker"),xlab="penalty",ylab="expected ransom")
#dev.off()

#plot.multiple(data.qre.r[,1],data.qre.r[,c(2,3)],c("target did not invest","target invested"),xlab="not-pay social norm",ylab="expected ransom")

#par.b.list<-list(c(par.b,0),c(par.b,1))
#data.qre.r<-calculate.multiple(qre.r.calc,par.b.list,key.index=8,range=c(0,50),n=20,par.setting=par.setting)

#plot.multiple(data.qre.r[,1],data.qre.r[,c(2,3)],c("target did not invest","target invested"),xlab="not-pay social norm",ylab="expected ransom")


# 
# qrr<-cbind(qrr1[,1],rep(0,nrow(qrr1)),rep(0,nrow(qrr1)))
# 
# for(i in 1:nrow(data.base)) {
#   if(data.base[i,"attack.decision"]>0) {
#     r<-max(1,data.base[i,"ransom"])
#     
#     old.payment.outcome<-ifelse(is.na(data.base[i,"old_payment_outcome"]),0,data.base[i,"old_payment_outcome"])
#     par.setting[7]<-old.payment.outcome
#     
#     ii<-ifelse(data.base[i,"attack.decision"]==1,data.base[i,"i_1"],data.base[i,"i_2"])
#     qrr.i<-qre.r(ii,par.setting=par.setting,par.behavior=par.b)
#     
#     qrr[r,3]<-qrr[r,3]+1
#     if(r<1) browser()
#     if(r>100) browser()
#     
#     qrr[,2]<-qrr[,2]+qrr.i[,2]
#   }
# }
# 
# qrr[,2]<-qrr[,2]/sum(qrr[,2])
# qrr[,3]<-qrr[,3]/sum(qrr[,3])
# 
# dimnames(qrr)[[2]]<-c("r","model","data")
# 
# qrr.plot<-plot.smooth(qrr,n=10)
# 
# ff<-paste("output/figure.",format(Sys.time(), "%Y_%m_%d_%H_%M_%S"),".",treat,".",model.list[model.index],".jpg",sep="")
# 
# m.lab<-paste(treat," ",model.list[model.index],sep="")
# plot.multiple(qrr.plot[,1],qrr.plot[,c(2,3)],c("model","data"),file=ff,main=m.lab,xlab="ransom",ylab="QRE prob")
# 

# att.summary<-cbind(c(0,1,2),c(0,0,0),c(0,0,0))
# for(i in 1:nrow(data.base)) {
#    inv.list<-c(data.base[i,"i_1"],data.base[i,"i_2"])
#    qra<-qre.att(inv.list,par.setting=par.setting,par.behavior=par.b)
#    att.summary[,2]<-att.summary[,2]+qra[,2]
#    att.summary[(att.summary[,1]==data.base[i,"attack.decision"]),3]<-att.summary[(att.summary[,1]==data.base[i,"attack.decision"]),3]+1
# }
#  
# att.summary[,2]<-att.summary[,2]/sum(att.summary[,2])
# att.summary[,3]<-att.summary[,3]/sum(att.summary[,3])
#  
# dimnames(att.summary)[[2]]<-c("att","model","data")
# 
# 

# qre.i.calc<-function(par,par.setting=NULL) {
# 
#   par.b<-par
#   
#   u.m1<-matrix(nrow=2,ncol=2)
#   u.m2<-matrix(nrow=2,ncol=2)
#   
#   for(i in 0:1) for(j in 0:1) {
#     u.m1[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
#     u.m2[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
#   }
#   
#   dimnames(u.m1)<-list(c("0","1"),c("0","1"))
#   dimnames(u.m2)<-list(c("0","1"),c("0","1"))
#   
#   qre.i<-find.qre(c(par.b[4],par.b[4]),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
#   # qre.i<-find.qre(c(0.01,0.01),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
#   
#   qi1<-qre.i[[1]]
#   qi2<-qre.i[[2]]
#   inv.summary<-c("no-no","no-yes","yes-no","yes-yes")
#   inv.summary<-cbind(inv.summary,c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2]))
#   
# #  browser()
#   
#   (qi1[2,2]+qi2[2,2])/2
# }
# 
# par.b.list<-list(c(par.b))
# data.qre.i.inv<-calculate.multiple(qre.i.calc,par.b.list,key.index=8,range=c(0,50),n=20,par.setting=par.setting)
# data.qre.i.pay<-calculate.multiple(qre.i.calc,par.b.list,key.index=6,range=c(0,50),n=20,par.setting=par.setting)
# 
# jpeg(file="output/sim.i.jpg")
# plot.multiple(data.qre.i.inv[,1],cbind(data.qre.i.inv[,2],data.qre.i.pay[,2]),c("invest util increase","not-pay util increase"),xlab="treatment utility",ylab="investment rate")
# dev.off()

qre.pay.calc<-function(par,par.behavior=NULL) {

  par.b<-par.behavior
  par.setting<-par

  u.m1<-matrix(nrow=2,ncol=2)
  u.m2<-matrix(nrow=2,ncol=2)

  for(i in 0:1) for(j in 0:1) {
    u.m1[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
    u.m2[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
  }

  dimnames(u.m1)<-list(c("0","1"),c("0","1"))
  dimnames(u.m2)<-list(c("0","1"),c("0","1"))

  qre.i<-find.qre(c(par.b[4],par.b[4]),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
  # qre.i<-find.qre(c(0.01,0.01),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")

  qi1<-qre.i[[1]]
  qi2<-qre.i[[2]]
  inv.summary<-c("no-no","no-yes","yes-no","yes-yes")
  inv.summary<-cbind(inv.summary,c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2]))
  inv.p<-c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2])


  qrr1<-qre.r(1,par.setting=par.setting,par.behavior=par.b)
  qrr0<-qre.r(0,par.setting=par.setting,par.behavior=par.b)

  qra00<-qre.att(c(0,0),par.setting=par.setting,par.behavior=par.b)
  qra01<-qre.att(c(0,1),par.setting=par.setting,par.behavior=par.b)
  qra10<-qre.att(c(1,0),par.setting=par.setting,par.behavior=par.b)
  qra11<-qre.att(c(1,1),par.setting=par.setting,par.behavior=par.b)

  prob.r0<-inv.p[1]*sum(qra00[2:3,2])+inv.p[2]*qra01[2,2]+inv.p[3]*qra10[3,2]
  prob.r1<-inv.p[2]*qra01[3,2]+inv.p[3]*qra10[2,2]+inv.p[4]*sum(qra11[2:3,2])

  prob.sum<-prob.r0+prob.r1

  prob.r0<-prob.r0/prob.sum
  prob.r1<-prob.r1/prob.sum

  pp1<-prob.pay(c(1:100),1,par.setting=par.setting,par.behavior=par.b)
  pp0<-prob.pay(c(1:100),0,par.setting=par.setting,par.behavior=par.b)

  pp<-sum(prob.r0*pp0*qrr0[,2]+prob.r1*pp1*qrr1[,2])
#  browser()
  pp
}

par.setting.list<-list(par.setting)
data.qre.pay.penalty<-calculate.multiple(qre.pay.calc,par.setting.list,key.index=6,range=c(0,100),n=20,par.behavior=par.b)
data.qre.pay.penalty<-cbind(data.qre.pay.penalty,rep(1,nrow(data.qre.pay.penalty)))
# data.qre.pay.pay<-calculate.multiple(qre.pay.calc,par.b.list,key.index=6,range=c(0,50),n=20,par.setting=par.setting)

#jpeg(file="output/sim.pay.jpg")
plot.multiple(data.qre.pay.penalty[,1],cbind(data.qre.pay.penalty[,2],data.qre.pay.penalty[,3]),
              c("boundedly rational defender","rational defender"),xlab="penalty",ylab="payment rate")
#dev.off()

# 
# qre.attack.calc<-function(par,par.setting=NULL) {
#   
#   par.b<-par
#   
#   u.m1<-matrix(nrow=2,ncol=2)
#   u.m2<-matrix(nrow=2,ncol=2)
#   
#   for(i in 0:1) for(j in 0:1) {
#     u.m1[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
#     u.m2[(i+1),(j+1)]<-u.inv(i,j,c(0,0),par.setting = par.setting,par.behavior = par.b)
#   }
#   
#   dimnames(u.m1)<-list(c("0","1"),c("0","1"))
#   dimnames(u.m2)<-list(c("0","1"),c("0","1"))
#   
#   qre.i<-find.qre(c(par.b[4],par.b[4]),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
#   # qre.i<-find.qre(c(0.01,0.01),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")
#   
#   qi1<-qre.i[[1]]
#   qi2<-qre.i[[2]]
#   inv.summary<-c("no-no","no-yes","yes-no","yes-yes")
#   inv.summary<-cbind(inv.summary,c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2]))
#   inv.p<-c(qi1[1,2]*qi2[1,2],qi1[1,2]*qi2[2,2],qi1[2,2]*qi2[1,2],qi1[2,2]*qi2[2,2])
#   
#   
#   qrr1<-qre.r(1,par.setting=par.setting,par.behavior=par.b)
#   qrr0<-qre.r(0,par.setting=par.setting,par.behavior=par.b)
#   
#   qra00<-qre.att(c(0,0),par.setting=par.setting,par.behavior=par.b)
#   qra01<-qre.att(c(0,1),par.setting=par.setting,par.behavior=par.b)
#   qra10<-qre.att(c(1,0),par.setting=par.setting,par.behavior=par.b)
#   qra11<-qre.att(c(1,1),par.setting=par.setting,par.behavior=par.b)
#   
#   prob.r0<-inv.p[1]*sum(qra00[2:3,2])+inv.p[2]*qra01[2,2]+inv.p[3]*qra10[3,2]
#   prob.r1<-inv.p[2]*qra01[3,2]+inv.p[3]*qra10[2,2]+inv.p[4]*sum(qra11[2:3,2])
#   
#   prob.sum<-prob.r0+prob.r1
#   
#   prob.r0<-prob.r0/prob.sum
#   prob.r1<-prob.r1/prob.sum
#   
#   pp1<-prob.pay(c(1:100),1,par.setting=par.setting,par.behavior=par.b)
#   pp0<-prob.pay(c(1:100),0,par.setting=par.setting,par.behavior=par.b)
#   
#   # pp<-sum(prob.r0*pp0*qrr0[,2]+prob.r1*pp1*qrr1[,2])
#   # browser()
#   
#   pp<-inv.p[1]*sum(qra00[2:3,2])+inv.p[2]*sum(qra01[2:3,2])+inv.p[3]*sum(qra10[2:3,2])+inv.p[4]*sum(qra11[2:3,2])
#   pp  
# } 
# 
# par.b.list<-list(c(par.b))
# data.qre.attack.inv<-calculate.multiple(qre.attack.calc,par.b.list,key.index=8,range=c(0,50),n=20,par.setting=par.setting)
# data.qre.attack.pay<-calculate.multiple(qre.attack.calc,par.b.list,key.index=6,range=c(0,50),n=20,par.setting=par.setting)
# 
# jpeg(file="output/sim.attack.jpg")
# plot.multiple(data.qre.attack.inv[,1],cbind(data.qre.attack.inv[,2],data.qre.attack.pay[,2]),c("invest util increase","not-pay util increase"),xlab="treatment utility",ylab="attack rate")
# dev.off()
# 

stopCluster(cluster)

time.end<-proc.time()
tt<-time.end[3]-time.start[3]
tt<-tt/60
cat("time elapsed = ",tt," minutes\n")
save.image()